#include "GridOperator.cuh"
#include "GridSolver.cuh"
#include "GridSampler.cuh"

// This kernel applies Dirichlet boundary condition on -X and +X planes
template<typename T>
__global__ void applyBoundaryXKernel(int nx, int ny, int nz, T bc, T *f)
{
	int y = blockIdx.x * blockDim.x + threadIdx.x;
	int z = blockIdx.y * blockDim.y + threadIdx.y;
	if(y < ny && z < nz)
	{
		int idx = y * nx + z * nx * ny;
		f[idx] = bc;
		f[idx + nx - 1] = bc;
	}
}

// This kernel applies Dirichlet boundary condition on -Y and +Y planes
template<typename T>
__global__ void applyBoundaryYKernel(int nx, int ny, int nz, T bc, T *f)
{
	int z = blockIdx.x * blockDim.x + threadIdx.x;
	int x = blockIdx.y * blockDim.y + threadIdx.y;
	if(z < nz && x < nx)
	{
		int idx = x + z * nx * ny;
		f[idx] = bc;
		f[idx + (ny - 1) * nx] = bc;
	}
}

// This kernel applies Dirichlet boundary condition on -Z and +Z planes
template<typename T>
__global__ void applyBoundaryZKernel(int nx, int ny, int nz, T bc, T *f)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	if(x < nx && y < ny)
	{
		int idx = x + y * nx;
		f[idx] = bc;
		f[idx + (nz - 1) * nx * ny] = bc;
	}
}

template<typename T>
void applyBoundaryX(int nx, int ny, int nz, T bc, T *f)
{
	dim3 gdim((ny - 1) / BLOCK_DIM + 1, (nz - 1) / BLOCK_DIM + 1);
	dim3 bdim(BLOCK_DIM, BLOCK_DIM, 1);
	applyBoundaryXKernel<T><<<gdim, bdim>>>(nx, ny, nz, bc, f);
	cutilCheckMsg("applyBoundaryXKernel() execution failed!");
}

template
void applyBoundaryX<float>(int nx, int ny, int nz, float bc, float *f);

template<typename T>
void applyBoundaryY(int nx, int ny, int nz, T bc, T *f)
{
	dim3 gdim((nz - 1) / BLOCK_DIM + 1, (nx - 1) / BLOCK_DIM + 1);
	dim3 bdim(BLOCK_DIM, BLOCK_DIM, 1);
	applyBoundaryYKernel<T><<<gdim, bdim>>>(nx, ny, nz, bc, f);
	cutilCheckMsg("applyBoundaryYKernel() execution failed!");
}

template
void applyBoundaryY<float>(int nx, int ny, int nz, float bc, float *f);

template<typename T>
void applyBoundaryZ(int nx, int ny, int nz, T bc, T *f)
{
	dim3 gdim((nx - 1) / BLOCK_DIM + 1, (ny - 1) / BLOCK_DIM + 1);
	dim3 bdim(BLOCK_DIM, BLOCK_DIM, 1);
	applyBoundaryZKernel<T><<<gdim, bdim>>>(nx, ny, nz, bc, f);
	cutilCheckMsg("applyBoundaryZKernel() execution failed!");
}

template
void applyBoundaryZ<float>(int nx, int ny, int nz, float bc, float *f);

// This kernel solves convection euqation on MAC grids
// The velocity values for convection are defined at face centers
// CX/CY/CZ describe how the scalar field to be convected is defined
// CX - if the scalar values are defined at cell centers along X-axis
// CY - if the scalar values are defined at cell centers along Y-axis
// CZ - if the scalar values are defined at cell centers along Z-axis
template<typename T, bool CX, bool CY, bool CZ>
__global__ void convectSemiLagrangeKernel(int nx, int ny, int nz, int z, T dx, T dt, T *u, T *v, T *w, T *f, T *fout)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	if(x < nx && y < ny)
	{
		T c = dt / dx;
		T x0 = x - sample<T>(CX ? nx + 1 : nx, CY ? ny : ny - 1, CZ ? nz : nz - 1, CX ? x + T(0.5) : x, CY ? y : y - T(0.5), CZ ? z : z - T(0.5), u) * c;
		T y0 = y - sample<T>(CX ? nx : nx - 1, CY ? ny + 1 : ny, CZ ? nz : nz - 1, CX ? x : x - T(0.5), CY ? y + T(0.5) : y, CZ ? z : z - T(0.5), v) * c;
		T z0 = z - sample<T>(CX ? nx : nx - 1, CY ? ny : ny - 1, CZ ? nz + 1 : nz, CX ? x : x - T(0.5), CY ? y : y - T(0.5), CZ ? z + T(0.5) : z, w) * c;
		int idx = x + y * nx + z * nx * ny;
		fout[idx] = sample<T>(nx, ny, nz, x0, y0, z0, f);
	}
}

template<typename T, bool CX, bool CY, bool CZ>
__global__ void convectMacCormackKernel(int nx, int ny, int nz, int z, T dx, T dt, T *u, T *v, T *w, T *f, T *fout)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	if(x < nx && y < ny)
	{
		T c = dt / dx;
		T x0 = x - sample<T>(CX ? nx + 1 : nx, CY ? ny : ny - 1, CZ ? nz : nz - 1, CX ? x + T(0.5) : x, CY ? y : y - T(0.5), CZ ? z : z - T(0.5), u) * c;
		T y0 = y - sample<T>(CX ? nx : nx - 1, CY ? ny + 1 : ny, CZ ? nz : nz - 1, CX ? x : x - T(0.5), CY ? y + T(0.5) : y, CZ ? z : z - T(0.5), v) * c;
		T z0 = z - sample<T>(CX ? nx : nx - 1, CY ? ny : ny - 1, CZ ? nz + 1 : nz, CX ? x : x - T(0.5), CY ? y : y - T(0.5), CZ ? z + T(0.5) : z, w) * c;
		T x1 = x0 + sample<T>(CX ? nx + 1 : nx, CY ? ny : ny - 1, CZ ? nz : nz - 1, CX ? x0 + T(0.5) : x0, CY ? y0 : y0 - T(0.5), CZ ? z0 : z0 - T(0.5), u) * c;
		T y1 = y0 + sample<T>(CX ? nx : nx - 1, CY ? ny + 1 : ny, CZ ? nz : nz - 1, CX ? x0 : x0 - T(0.5), CY ? y0 + T(0.5) : y0, CZ ? z0 : z0 - T(0.5), v) * c;
		T z1 = z0 + sample<T>(CX ? nx : nx - 1, CY ? ny : ny - 1, CZ ? nz + 1 : nz, CX ? x0 : x0 - T(0.5), CY ? y0 : y0 - T(0.5), CZ ? z0 + T(0.5) : z0, w) * c;
		T x2 = x0 + (x - x1) * T(0.5);
		T y2 = y0 + (y - y1) * T(0.5);
		T z2 = z0 + (z - z1) * T(0.5);
		int idx = x + y * nx + z * nx * ny;
		fout[idx] = sample<T>(nx, ny, nz, x2, y2, z2, f);
	}
}

template<typename T>
void convect(int nx, int ny, int nz, T dx, T dt, T *u, T *v, T *w, T *f, T *fout, const FieldType type, const NumericalScheme scheme)
{
	dim3 gdim((nx - 1) / BLOCK_DIM + 1, (ny - 1) / BLOCK_DIM + 1);
	dim3 bdim(BLOCK_DIM, BLOCK_DIM, 1);
	switch(scheme)
	{
	case SemiLagrange:
		switch(type)
		{
		case CellCenter:
			for(int z = 0; z < nz; ++z)
				convectSemiLagrangeKernel<T, true, true, true><<<gdim, bdim>>>(nx, ny, nz, z, dx, dt, u, v, w, f, fout);
			break;
		case FaceCenterX:
			for(int z = 0; z < nz; ++z)
				convectSemiLagrangeKernel<T, false, true, true><<<gdim, bdim>>>(nx, ny, nz, z, dx, dt, u, v, w, f, fout);
			break;
		case FaceCenterY:
			for(int z = 0; z < nz; ++z)
				convectSemiLagrangeKernel<T, true, false, true><<<gdim, bdim>>>(nx, ny, nz, z, dx, dt, u, v, w, f, fout);
			break;
		case FaceCenterZ:
			for(int z = 0; z < nz; ++z)
				convectSemiLagrangeKernel<T, true, true, false><<<gdim, bdim>>>(nx, ny, nz, z, dx, dt, u, v, w, f, fout);
			break;
		}
		cutilCheckMsg("convectSemiLagrangeKernel() execution failed!");
		break;
	case MacCormack:
		switch(type)
		{
		case CellCenter:
			for(int z = 0; z < nz; ++z)
				convectMacCormackKernel<T, true, true, true><<<gdim, bdim>>>(nx, ny, nz, z, dx, dt, u, v, w, f, fout);
			break;
		case FaceCenterX:
			for(int z = 0; z < nz; ++z)
				convectMacCormackKernel<T, false, true, true><<<gdim, bdim>>>(nx, ny, nz, z, dx, dt, u, v, w, f, fout);
			break;
		case FaceCenterY:
			for(int z = 0; z < nz; ++z)
				convectMacCormackKernel<T, true, false, true><<<gdim, bdim>>>(nx, ny, nz, z, dx, dt, u, v, w, f, fout);
			break;
		case FaceCenterZ:
			for(int z = 0; z < nz; ++z)
				convectMacCormackKernel<T, true, true, false><<<gdim, bdim>>>(nx, ny, nz, z, dx, dt, u, v, w, f, fout);
			break;
		}
		cutilCheckMsg("convectMacCormackKernel() execution failed!");
		break;
	}
}

template
void convect<float>(int nx, int ny, int nz, float dx, float dt, float *u, float *v, float *w, float *f, float *fout, const FieldType type, const NumericalScheme scheme);

// This kernel is designed for solving Poisson/heat euqations with zero Neumann boundary conditions
// These equations has a general form: (aI + bL) u = f where L is the Laplacian operator (SPD)
// Red-Black Gauss-Seidel iteration is supported through template arguments
// IR - ignore red points
// IB - ignore black points
// a - use 0 for Poisson equation and 1 for heat equation
// b - use 1 for Poisson equation and k\Delta t for heat equation
// f - right hand side vector
// u - current unknown vector before iteration
// uout - new unknown vector after iteration
template<typename T, bool IR, bool IB>
__global__ void iterateKernel(int nx, int ny, int nz, int z, T dx, T a, T b, T *f, T *u, T *uout)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	if(IR) { if((x + y + z) & 1 == 0) return; }
	if(IB) { if((x + y + z) & 1 == 1) return; }
	if(x < nx && y < ny)
	{
		int idx = x + y * nx + z * nx * ny;
		T aij = b / (dx * dx);
		T aii(a), uij(0);
		if(x > 0)
		{ aii += aij; uij += u[idx - 1]; }
		if(x < nx - 1)
		{ aii += aij; uij += u[idx + 1]; }
		if(y > 0)
		{ aii += aij; uij += u[idx - nx]; }
		if(y < ny - 1)
		{ aii += aij; uij += u[idx + nx]; }
		if(z > 0)
		{ aii += aij; uij += u[idx - nx * ny]; }
		if(z < nz - 1)
		{ aii += aij; uij += u[idx + nx * ny]; }
		uout[idx] = (f[idx] + uij * aij) / aii;
	}
}

template<typename T>
void solveHeatEquation(int nx, int ny, int nz, T dx, T dt, T kappa, T *u, T *uout, int its)
{
	dim3 gdim((nx - 1) / BLOCK_DIM + 1, (ny - 1) / BLOCK_DIM + 1);
	dim3 bdim(BLOCK_DIM, BLOCK_DIM, 1);
	CUDA_SAFE_CALL(cudaMemcpy(uout, u, nx * ny * nz * sizeof(T), cudaMemcpyDeviceToDevice));
	T kdt = kappa * dt;
	for(int it = 0; it < its; ++it)
	{
		for(int z = 0; z < nz; ++z)
			iterateKernel<T, true, false><<<gdim, bdim>>>(nx, ny, nz, z, dx, T(1), kdt, u, uout, uout);
		cutilCheckMsg("iterateKernel() execution failed!");

		for(int z = 0; z < nz; ++z)
			iterateKernel<T, false, true><<<gdim, bdim>>>(nx, ny, nz, z, dx, T(1), kdt, u, uout, uout);
		cutilCheckMsg("iterateKernel() execution failed!");
	}
}

template
void solveHeatEquation<float>(int nx, int ny, int nz, float dx, float dt, float kappa, float *u, float *uout, int its);

template<typename T>
void solvePoissonEquation(int nx, int ny, int nz, T dx, T *f, T *u, int its)
{
	dim3 gdim((nx - 1) / BLOCK_DIM + 1, (ny - 1) / BLOCK_DIM + 1);
	dim3 bdim(BLOCK_DIM, BLOCK_DIM, 1);
	CUDA_SAFE_CALL(cudaMemset(u, 0, nx * ny * nz * sizeof(T)));
	for(int it = 0; it < its; ++it)
	{
		for(int z = 0; z < nz; ++z)
			iterateKernel<T, true, false><<<gdim, bdim>>>(nx, ny, nz, z, dx, T(0), T(1), f, u, u);
		cutilCheckMsg("iterateKernel() execution failed!");

		for(int z = 0; z < nz; ++z)
			iterateKernel<T, false, true><<<gdim, bdim>>>(nx, ny, nz, z, dx, T(0), T(1), f, u, u);
		cutilCheckMsg("iterateKernel() execution failed!");
	}
}

template
void solvePoissonEquation<float>(int nx, int ny, int nz, float dx, float *f, float *u, int its);

// This kernel performs velocity projection using pressure gradient on MAC grids
// Here (nx, ny, nz) is the cell dimension
template<typename T>
__global__ void projectKernel(int nx, int ny, int nz, int z, T dx, T *u, T *v, T *w, T *p)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	if(x < nx && y < ny)
	{
		int idx = x + y * nx + z * nx * ny;
		T rdx = T(1) / dx;
		T p0 = p[idx];
		if(x > 0) u[idx + y + z * ny] -= (p0 - p[idx - 1]) * rdx;
		if(y > 0) v[idx + z * nx] -= (p0 - p[idx - nx]) * rdx;
		if(z > 0) w[idx] -= (p0 - p[idx - nx * ny]) * rdx;
	}
}

template<typename T>
void project(int nx, int ny, int nz, T dx, T *u, T *v, T *w, T *p)
{
	dim3 gdim((nx - 1) / BLOCK_DIM + 1, (ny - 1) / BLOCK_DIM + 1);
	dim3 bdim(BLOCK_DIM, BLOCK_DIM, 1);
	for(int z = 0; z < nz; ++z)
		projectKernel<T><<<gdim, bdim>>>(nx, ny, nz, z, dx, u, v, w, p);
	cutilCheckMsg("projectKernel() execution failed!");
}

template
void project<float>(int nx, int ny, int nz, float dx, float *u, float *v, float *w, float *p);

// This kernel apply gravity force to y-component of velocity fields on MAC grids
// v += g * (f - famb) * dt
template<typename T>
__global__ void gravityKernel(int nx, int ny, int nz, int z, T dx, T dt, T g, T famb, T *f, T *v)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	if(x < nx && y < ny)
	{
		int idx = x + y * nx + z * nx * ny;
		int idx1 = x + min(max((y - 1), 0), ny - 2) * nx + z * nx * (ny - 1);
		int idx2 = x + min(max(y, 0), ny - 2) * nx + z * nx * (ny - 1);
		v[idx] += ((f[idx1] + f[idx2]) * T(0.5) - famb) * g * dt;
	}
}

template<typename T>
void gravity(int nx, int ny, int nz, T dx, T dt, T g, T famb, T *f, T *v)
{
	dim3 gdim((nx - 1) / BLOCK_DIM + 1, (ny - 1) / BLOCK_DIM + 1);
	dim3 bdim(BLOCK_DIM, BLOCK_DIM, 1);
	for(int z = 0; z < nz; ++z)
		gravityKernel<T><<<gdim, bdim>>>(nx, ny, nz, z, dx, dt, g, famb, f, v);
	cutilCheckMsg("gravityKernel() execution failed!");
}

template
void gravity(int nx, int ny, int nz, float dx, float dt, float g, float famb, float *f, float *v);

// This kernel evaluates the vorticity confinement fields on MAC grids
// The values of vector fields are defined at cell centers
template<typename T>
__global__ void turbulenceKernel(int nx, int ny, int nz, int z, T dx, T eps, T *cx, T *cy, T *cz, T *cn, T *fx, T *fy, T *fz)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	if(x < nx && y < ny)
	{
		int nxy = nx * ny;
		int idx = x + y * nx + z * nxy;
		int id1 = (x == 0 ? idx : idx - 1);
		int id2 = (x == nx - 1 ? idx : idx + 1);
		int id3 = (y == 0 ? idx : idx - nx);
		int id4 = (y == ny - 1 ? idx : idx + nx);
		int id5 = (z == 0 ? idx : idx - nxy);
		int id6 = (z == nz - 1 ? idx : idx + nxy);
		T nxi = cn[id2] - cn[id1];
		T nyi = cn[id4] - cn[id3];
		T nzi = cn[id6] - cn[id5];
		T rni = T(1) / max(sqrt(nxi * nxi + nyi * nyi + nzi * nzi), T(1e-6));
		nxi *= rni;
		nyi *= rni;
		nzi *= rni;
		T cxi = cx[idx];
		T cyi = cy[idx];
		T czi = cz[idx];
		T edx = eps * dx;
		fx[idx] = (nyi * czi - nzi * cyi) * edx;
		fy[idx] = (nzi * cxi - nxi * czi) * edx;
		fz[idx] = (nxi * cyi - nyi * cxi) * edx;
	}
}

template<typename T>
void turbulence(int nx, int ny, int nz, T dx, T eps, T *cx, T *cy, T *cz, T *cn, T *fx, T *fy, T *fz)
{
	dim3 gdim((nx - 1) / BLOCK_DIM + 1, (ny - 1) / BLOCK_DIM + 1);
	dim3 bdim(BLOCK_DIM, BLOCK_DIM, 1);
	for(int z = 0; z < nz; ++z)
		turbulenceKernel<T><<<gdim, bdim>>>(nx, ny, nz, z, dx, eps, cx, cy, cz, cn, fx, fy, fz);
	cutilCheckMsg("turbulenceKernel() execution failed!");
}

template
void turbulence<float>(int nx, int ny, int nz, float dx, float eps, float *cx, float *cy, float *cz, float *cn, float *fx, float *fy, float *fz);

// This kernel accelerates velocity fields from cell centered forces on MAC grids
// The values of velocity fields are defined at face centers
template<typename T>
__global__ void accelerateKernel(int nx, int ny, int nz, int z, T dx, T dt, T *fx, T *fy, T *fz, T *u, T *v, T *w)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	if(x < nx && y < ny)
	{
		int nxy = nx * ny;
		int idx = x + y * nx + z * nxy;
		int idxu = idx + y + z * ny;
		int idxv = idx + z * nx;
		int idxw = idx;
		T hdt = dt * T(0.5);
		T fxdt = fx[idx] * hdt;
		T fydt = fy[idx] * hdt;
		T fzdt = fz[idx] * hdt;
		u[idxu] += fxdt;
		u[idxu + 1] += fxdt;
		v[idxv] += fydt;
		v[idxv + nx] += fydt;
		w[idxw] += fzdt;
		w[idxw + nxy] += fzdt;
	}
}

template<typename T>
void accelerate(int nx, int ny, int nz, T dx, T dt, T *fx, T *fy, T *fz, T *u, T *v, T *w)
{
	dim3 gdim((nx - 1) / BLOCK_DIM + 1, (ny - 1) / BLOCK_DIM + 1);
	dim3 bdim(BLOCK_DIM, BLOCK_DIM, 1);
	for(int z = 0; z < nz; ++z)
		accelerateKernel<T><<<gdim, bdim>>>(nx, ny, nz, z, dx, dt, fx, fy, fz, u, v, w);
	cutilCheckMsg("accelerateKernel() execution failed!");
}

template
void accelerate<float>(int nx, int ny, int nz, float dx, float dt, float *fx, float *fy, float *fz, float *u, float *v, float *w);

// This kernel solves the reaction equation on MAC grids
// The values of fields are defined at cell centers
// c - cooling coefficient
// f - fuel
// s - soot
// t - temperature
template<typename T>
__global__ void reactionKernel(int nx, int ny, int nz, int z, T dx, T dt, T burn, T cool, T tamb, T tmax, T *f, T *s, T *t)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	if(x < nx && y < ny)
	{
		int idx = x + y * nx + z * nx * ny;
		T fi = f[idx];
		T si = s[idx];
		T ti = t[idx] * T(0.99);
		int power = 1;
		if(x > 0 && f[idx - 1] < T(0.001)) ++power;
		if(x < nx - 1 && f[idx + 1] < T(0.001)) ++power;
		if(y > 0 && f[idx - nx] < T(0.001)) ++power;
		if(y < ny - 1 && f[idx + nx] < T(0.001)) ++power;
		if(z > 0 && f[idx - nx * ny] < T(0.001)) ++power;
		if(z < nz - 1 && f[idx + nx * ny] < T(0.001)) ++power;
		burn *= T(power) / T(7);
		T transfer = min(burn * dt, fi);
		f[idx] = fi - transfer;
		s[idx] = si + transfer;
		T c1 = ti - tamb;
		c1 *= c1 * c1;
		T c2 = tmax - tamb;
		c2 *= c2;
		c2 *= c2;
		t[idx] = (tamb + pow(T(1) / c1 + T(3) * cool * dt / c2, -T(1) / T(3)) + transfer) / (T(1) + T(1) * dt);
	}
}

template<typename T>
void reaction(int nx, int ny, int nz, T dx, T dt, T burn, T cool, T tamb, T tmax, T *f, T *s, T *t)
{
	dim3 gdim((nx - 1) / BLOCK_DIM + 1, (ny - 1) / BLOCK_DIM + 1);
	dim3 bdim(BLOCK_DIM, BLOCK_DIM, 1);
	for(int z = 0; z < nz; ++z)
		reactionKernel<T><<<gdim, bdim>>>(nx, ny, nz, z, dx, dt, burn, cool, tamb, tmax, f, s, t);
	cutilCheckMsg("reactionKernel() execution failed!");
}

template
void reaction<float>(int nx, int ny, int nz, float dx, float dt, float burn, float cool, float tamb, float tmax, float *f, float *s, float *t);

template<typename T>
__global__ void convectParticlesKernel(int nx, int ny, int nz, T dx, T dt, T *u, T *v, T *w, int np, float3 *p)
{
	int idx = blockIdx.x * blockDim.x + threadIdx.x;
	if(idx < np)
	{
		float3 pi = p[idx];
		T rdx = T(1) / dx;
		T xi = pi.x * rdx;
		T yi = pi.y * rdx;
		T zi = pi.z * rdx;
		T ui = sample(nx + 1, ny, nz, xi, yi - T(0.5), zi - T(0.5), u);
		T vi = sample(nx, ny + 1, nz, xi - T(0.5), yi, zi - T(0.5), v);
		T wi = sample(nx, ny, nz + 1, xi - T(0.5), yi - T(0.5), zi, w);
		pi.x += ui * dt;
		pi.y += vi * dt;
		pi.z += wi * dt;
		p[idx] = pi;
	}
}

template<typename T>
void convectParticles(int nx, int ny, int nz, T dx, T dt, T *u, T *v, T *w, int np, float3 *p)
{
	const int bdim2 = BLOCK_DIM * BLOCK_DIM;
	dim3 gdim((np - 1) / bdim2 + 1, 1);
	dim3 bdim(bdim2, 1, 1);
	convectParticlesKernel<T><<<gdim, bdim>>>(nx, ny, nz, dx, dt, u, v, w, np, p);
	cutilCheckMsg("convectParticlesKernel() execution failed!");
}

template
void convectParticles(int nx, int ny, int nz, float dx, float dt, float *u, float *v, float *w, int np, float3 *p);
